\section{用SageMath计算举例}

\begin{frame}[fragile]{原根}

  \mintinline{sage}{primitive_root(n)} 提供了模 $n=2, 4, p^s$ 或 $2p^s$ 的一个原根；
  且如果$n$为素数，这给出了最小原根。
我们可以写个函数来找到所有的模$m$原根。

\begin{minted}[texcl]{sage}
sage: primitive_root(43)
3
sage: primitive_root(43^2)
3
sage: primitive_root(2*43^2)
3
sage: def all_primitive_roots(n):  # 求模 $n$ 的所有原根的函数
....:     a_primitive_root = primitive_root(n)
....:     phi_n = euler_phi(n)
....:     R = IntegerModRing(n)
....:     all = []
....:     for i in range(1, phi_n):
....:         if Integer(i).gcd(phi_n) == 1:
....:             all.append(R(a_primitive_root)^i)
....:     print(all)
....: 
sage: all_primitive_roots(19)
[2, 13, 14, 15, 3, 10]
\end{minted}
\end{frame}

\begin{frame}[fragile]{$\ZZ/m\ZZ$的单位群$U_m$的结构、元素的阶}
  要知道$U_m$的结构，可调用环的 \mintinline{sage}{unit_group()} 方法。
  求$a\in U_m$的阶，可用 \mintinline{sage}{multiplicative_order()} 方法。

\begin{minted}[texcl]{sage}
sage: R=IntegerModRing(43)
sage: R.unit_group()  # $U_{43}\cong C_{42}$
Multiplicative Abelian group isomorphic to C42
sage: R=IntegerModRing(72)  
sage: R.unit_group()  # $U_{72}\cong C_2\times C_2\times C_6$
Multiplicative Abelian group isomorphic to C2 x C2 x C6
sage: R=IntegerModRing(2^3*3^4*5^2)  
sage: R.unit_group()  # $U_{2^3\times 3^4\times 5^2}\cong C_2\times C_2\times C_{54}\times C_{20}$
Multiplicative Abelian group isomorphic to C2 x C2 x C54 x C20   
sage: R=IntegerModRing(50)
sage: R.unit_group()  # $U_{50}\cong C_{20}$
Multiplicative Abelian group isomorphic to C20
sage: R(13).multiplicative_order()  # $\ord_{50}13=20$
20
sage: R(9).multiplicative_order()  # $\ord_{50}9=10$
10
\end{minted}
  \end{frame}

\begin{frame}[fragile]{离散对数}

  求解离散对数问题，即找$x$使得 $b^x\equiv a\left( \mod m \right)$, 可以用 \mintinline{sage}{log} 方法。
\begin{minted}[texcl]{sage}
sage: m = 2*5^2
sage: R = IntegerModRing (m)  # 定义环$\ZZ/50\ZZ$
sage: g = primitive_root (m)  # 一个原根 $g=27$ 
sage: g
27
sage: R(19).log(g)  # 计算$19\left( \mod 50 \right)$相对于原根$g$的指标 $\ind_g(19)$
18
sage: R(g)^18  # 验证指标确为 $18$
19
sage: R(19).log(29)  # 求 $x$ 使得 $29^x\equiv 19\left( \mod 50 \right)$
9
\end{minted}

\end{frame}

\begin{frame}[fragile]{$n$次剩余问题}
  \mintinline{sage}{is_square()} 可用于检验是否为二次剩余。
  \mintinline{sage}{nth_root(n)}方法可用于找$n$次根 (只要有)。

\begin{minted}[texcl]{sage}
sage: R = IntegerModRing(11)  # 定义环 $\ZZ/11\ZZ$
sage: R(3).is_square()  # $\bar{3}\in U_{11}$是否为二次剩余
True
sage: def all_squares(n):  # 求出模 $n$ 的所有二次剩余
....:     R = IntegerModRing(n)
....:     all = []
....:     for i in range(1, n):
....:         if Integer(i).gcd(n) == 1 and R(i).is_square():
....:             all.append(i)
....:     print(all)
....: 
sage: all_squares(23)  # 所有的模$23$二次剩余
[1, 2, 3, 4, 6, 8, 9, 12, 13, 16, 18]
sage: R=IntegerModRing(101)
sage: R(23).nth_root(2)
86
sage: R(23).nth_root(3)
45
sage: R(23).nth_root(5)
# 这里省略了错误信息
ValueError: no nth root
\end{minted}
\end{frame}

\begin{frame}[fragile]

  我们可以计算出所有$U_m$中元素的$n$次幂来得到所有模$m$的$n$次剩余。
\begin{minted}[texcl]{sage}
sage: def all_n_th_power(m, n):  # 求所有的模$m$的$n$次剩余
....:     R = IntegerModRing(m)
....:     all = []
....:     for i in range(1, m):
....:         if Integer(i).gcd(m) == 1:
....:             all.append(R(i)^n)
....:     print(set(all))
....: 
sage: all_n_th_power(11, 2)  # 模$11$的$2$次剩余
{1, 3, 4, 5, 9}
sage: all_n_th_power(11, 3)  # 模$11$的$3$次剩余
{1, 2, 3, 4, 5, 6, 7, 8, 9, 10}
sage: all_n_th_power(11, 4)  # 模$11$的$4$次剩余
{1, 3, 4, 5, 9}
sage: all_n_th_power(11, 5)  # 模$11$的$5$次剩余
{1, 10}
sage: all_n_th_power(12, 2)  # 模$12$的$2$次剩余
{1}
sage: all_n_th_power(12, 5)  # 模$12$的$5$次剩余
{1, 11, 5, 7}
\end{minted}
\end{frame}


\begin{frame}[fragile]{高次同余方程求解}
\begin{minted}[texcl]{sage}
sage: R.<x>=IntegerModRing(11)[]  # 定义环 $\ZZ/11\ZZ[x]$
sage: f=x^2-7
sage: f.roots()  # $x^2-7\in \ZZ/11\ZZ[x]$ 无根
[]
sage: R(7).is_square()  
False
sage: f=x^3-3
sage: f.roots()  # $x^3-3\in \ZZ/11\ZZ[x]$ 无根
[]
sage: f=x^3-2*x^2-3  
sage: f.roots()  # $x^3-2x-3\in \ZZ/11\ZZ[x]$ 只有一个$1$重根$7$
[(7, 1)]
sage: R.<x>=IntegerModRing(25)[]
sage: f=x^2-6
sage: f.roots(multiplicities=False)  # $x^2-6\in \ZZ/25\ZZ[x]$有两个根$9,16$
[9, 16]
sage: f=x^3-7
sage: f.roots(multiplicities=False)  # $x^3-7\in \ZZ/25\ZZ[x]$有一个根$18$
[18]
\end{minted}
\end{frame}

\begin{frame}[fragile]{整系数多项式方程的$p$-进制解}
  \mintinline{sage}{Zp} 类用于构造$p$-进制环$\ZZ_p=\varprojlim \ZZ/p^n\ZZ$.
  $p$-进制整数可想成级数$\sum_{i=0}^\infty a_i p^i$, 其中$0\leqslant a_i <p$.  
  在SageMath中可以指定精度计算到你需要的地步。
  若精度为$n$, 则$p$-进制环中元素会表示为
  \[
    a_0+a_1 p + \cdots + a_{n-1} p^{n-1} + O(p^{n}).
  \]
  关于SageMath中$p$-进制环的介绍，参见\href{https://doc.sagemath.org/html/en/reference/padics/sage/rings/padics/tutorial.html}{这里}。

\begin{minted}[texcl]{sage}
sage: # 如下定义精度为$10$的$7$-进制环 $BR$
sage: BR = Zp(7, prec = 10, type = 'fixed-mod', print_mode = 'series')  
sage: R.<x> = BR[]  # 定义多项式环 $R=BR[x]$
sage: f=x^2-2
sage: f.roots()  # 求 $x^2-2$ 的$7$-进制解（精度为$10$）
[(3 + 7 + 2*7^2 + 6*7^3 + 7^4 + 2*7^5 + 7^6 + 2*7^7 + 4*7^8 + 6*7^9, 1),
 (4 + 5*7 + 4*7^2 + 5*7^4 + 4*7^5 + 5*7^6 + 4*7^7 + 2*7^8, 1)]
\end{minted}

应用Hensel引理的证明中得到的公式
\[
  x_{s} \equiv x_{s-1}-f\left(x_{s-1}\right) / f^{\prime}\left(x_{s-1}\right)\left( \mod p^{s+1} \right)
\]
我们容易自己实现一个提升的函数。
可把下面的代码放入文本文件 (如 \mintinline{sage}{hensel_lift.sage})中，
然后通过 SageMath 中的 \mintinline{sage}{load} 函数加载。
\end{frame}

\begin{frame}[fragile]

\begin{minted}[texcl]{sage}
def my_hensel_lift(p, f, x0, prec=6):
    x0 = x0 % p
    if f.subs(x=x0) % p != 0:
        raise Exception('x0 不是所给多项式 f 模 p 的根')
    f_diff = f.differentiate()
    if f_diff.subs(x=x0) % p == 0:
        # 此时，参见 Rosen《初等数论及其应用》中定理 4.15
        raise Exception('f\' 在 x0 处模 p 为 0, 此函数未考虑')   
    R = IntegerModRing(p)
    coslope = Integer(R(f_diff.subs(x=x0))^(-1))
    coefficients, p_adic_solution = [x0], [x0]
    for s in range(1, prec):
        step = -f.subs(x=p_adic_solution[-1]) * coslope % p^(s+1)
        p_adic_solution.append(p_adic_solution[-1]+step)
        a_s = step / p^s
        coefficients.append(a_s)
    last_x_s = str(x0)
    for s in range(1, prec):
        x_s = ' + '.join([last_x_s, f'{coefficients[s]}*{p}^{s}'])
        last_x_s = x_s
        print(f'x{s} = {x_s}')
\end{minted}

\end{frame}

\begin{frame}[fragile]
  如下加载并调用：
\begin{minted}[texcl]{sage}
sage: load('hensel_lift.sage')
sage: Zx.<x>=ZZ[]
sage: my_hensel_lift(7, x^2-2, 3, 9)
x1 = 3 + 1*7^1
x2 = 3 + 1*7^1 + 2*7^2
x3 = 3 + 1*7^1 + 2*7^2 + 6*7^3
x4 = 3 + 1*7^1 + 2*7^2 + 6*7^3 + 1*7^4
x5 = 3 + 1*7^1 + 2*7^2 + 6*7^3 + 1*7^4 + 2*7^5
x6 = 3 + 1*7^1 + 2*7^2 + 6*7^3 + 1*7^4 + 2*7^5 + 1*7^6
x7 = 3 + 1*7^1 + 2*7^2 + 6*7^3 + 1*7^4 + 2*7^5 + 1*7^6 + 2*7^7
x8 = 3 + 1*7^1 + 2*7^2 + 6*7^3 + 1*7^4 + 2*7^5 + 1*7^6 + 2*7^7 + 4*7^8
sage: my_hensel_lift(7, x^2-2, 4, 9)
x1 = 4 + 5*7^1
x2 = 4 + 5*7^1 + 4*7^2
x3 = 4 + 5*7^1 + 4*7^2 + 0*7^3
x4 = 4 + 5*7^1 + 4*7^2 + 0*7^3 + 5*7^4
x5 = 4 + 5*7^1 + 4*7^2 + 0*7^3 + 5*7^4 + 4*7^5
x6 = 4 + 5*7^1 + 4*7^2 + 0*7^3 + 5*7^4 + 4*7^5 + 5*7^6
x7 = 4 + 5*7^1 + 4*7^2 + 0*7^3 + 5*7^4 + 4*7^5 + 5*7^6 + 4*7^7
x8 = 4 + 5*7^1 + 4*7^2 + 0*7^3 + 5*7^4 + 4*7^5 + 5*7^6 + 4*7^7 + 2*7^8
\end{minted}
\end{frame}

\begin{frame}[fragile]{算术函数}
\mintinline{sage}{moebius(n)} 用于计算M\"obius函数 $\mu(n)$. 

\begin{minted}[texcl]{sage}
sage: moebius(1)
1
sage: moebius(12); moebius(81)
0
0
sage: moebius(13); moebius(3*5*7)
-1
-1
sage: moebius(2*3*5*7)
1
\end{minted}

$\sigma_k(n)$按定义为正整数$n$的所有正整数因子的$k$次幂之和，其封闭的表达式为
\[
  \sigma_k(p_1^{e_1} \cdots p_s^{e_s})=\prod_{i=1}^s (1+ p_i^k+ p_i^{2k}\cdots +p_i^{e_ik}).
\]
\mintinline{sage}{sigma(n,k)} 用于计算$\sigma_k(n)$：
\begin{minted}[texcl]{sage}
sage: sigma(100,0);sigma(100,2)
9
13671
\end{minted}


\end{frame}
